function [outputresultmap outputresulterrormap] = dGemricQmapkneeregister(basedir,serie,slice)
% ratheartregister: run registration between multiple days
%
% Usage:  outputresultmap = ratheartregister(fixedstudydir, movingstudydir)
%
% Where: 
%   fixedstudydir, movingstudydir are the T2 directories of day 1, and
% a next day, respectively.
%   outputresultmap = the name of the registered (moving) R2map dcm.
%

% Load all images in directory and saves as one 3D nifti image
% In the basedir should be all the volunteerdirs, 'locations.mat' and
% 'Estimated_parameters.mat'
% Henk Smit

% Load image locations
cd(basedir)
load('locations.mat')

% Fixed image dir
for i=1:size(locations(serie).loc,1)-1
    QmapsFixed(i).dir = locations(serie).loc(i+1,slice).dir;
    QmapsFixed(i).file = locations(serie).loc(i+1,slice).file;
    QmapsFixed(i).slidenr = str2num(locations(serie).loc(i+1,slice).file(1:4));
    fixedims(i).im = fullfile(QmapsFixed(i).dir, QmapsFixed(i).file);
end
cd(QmapsFixed(i).dir);
cd ..

% Volunteer serie dir
volunteerdir=pwd;
d=dir;

% Number of Qmap inversion times
nrimages = 5;

% Postprocessing directory
postprocdir = fullfile(volunteerdir,'Postprocessing');
if ~exist(postprocdir, 'dir')   
      [success, message, messageid] = mkdir(postprocdir);
      if ~success    
         error( 'ratheartregister:invalidoutputdir', message );     
      end 
end

%Write dicom images in 1 3D nifti image
% for i=3:size(d,1)
%     if size(char(regexp(d(i).name,'SPGR', 'match')),2)>=1
%         directory = d(i).name;
%         [tstimg, info_structure] = loaddata([directory filesep], 1 , 1.5);
%         tstimg = squeeze(tstimg);
%         T = affineTransformFromDicominfo( info_structure );
%         [nii] = make_niiT(uint16(tstimg), T);
%         save_nii( nii, [directory, filesep,'3D.nii'])
%     end
% end

% Moving images

for j=1:size(QmapsFixed,2)   
    for i=1:nrimages-1
        filenr=QmapsFixed(j).slidenr+i;
        if filenr>9       
            movingims(i,j).name = [QmapsFixed(j).dir,'00',num2str(filenr),'.dcm'];
        else
            movingims(i,j).name = [QmapsFixed(j).dir,'000',num2str(filenr),'.dcm'];
        end
    end
end

% find elastix.exe
currentdir= fileparts( which('dGEmricIRregister.m') );
cd(currentdir)
cd ..
elastixdir = [pwd filesep 'elastix'];



% Parameter file
par = fullfile(elastixdir, 'parameters_RigidLMI2D.txt');

for k=1:size(movingims,2)
        % Registration mask
        tstinfo = dicominfo(movingims(1,k).name);
        tstimg = dicomread(tstinfo);
        
        imagemask = zeros( size(tstimg) );
        xmask = round(0.16*size(tstimg,1)):round(0.83*size(tstimg,1)); %square that covers the bone
        ymask = round(0.117*size(tstimg,1)):round(0.83*size(tstimg,1));
        imagemask(xmask,ymask,:) = 1;

        % Write mask to outputfolder
        regmaskdir = fullfile(postprocdir, 'regmaskqmap.dcm');
        dicomwrite(imagemask,regmaskdir,tstinfo);
        fixedmask = regmaskdir;
        movingmask = regmaskdir;
        
        for r=1:size(movingims,1) % Perform registration for all moving images
            display(['Performing registration ',num2str(r),'/',num2str(size(movingims,1)),' of protocol ',num2str(k),'/',num2str(size(movingims,2))])

            % Output directory
            outdir = fullfile(volunteerdir,'Postprocessing','RegistrationsQMAP',movingims(r,k).name(end-16:end-9));
            if ~exist(outdir, 'dir')   
                  [success, message, messageid] = mkdir(outdir);
                  if ~success    
                     error( 'ratheartregister:invalidoutputdir', message );     
                  end 
            end
            
            fixedim = fixedims(1,k).im;
            movingim = movingims(r,k).name;

            cmdline1 = ['"' fullfile(elastixdir, 'elastix') ...
            '" -f "' fixedim '" -m "' movingim ...
            '" -fMask "' fixedmask '" -mMask "' movingmask ...
            '" -p "' par '" -out "' outdir '"'];

            % Call Elastix and perform registration
            [status1,result1] = system(cmdline1);
            
            %copy fixed image to outdir
            if r==1
                outfixedim = fullfile(outdir,fixedims(1,k).im(end-7:end));
                copyfile(fixedim,outfixedim);
            end
  
            %rename result to original name
            outputfile = fullfile(outdir,'result.0.dcm');
            resultfile = fullfile(outdir,movingims(r,k).name(end-7:end));
            movinginfo = dicominfo(movingim); 
            resultdicom = dicomread(outputfile);
            dicomwrite(resultdicom,resultfile,movinginfo,'CreateMode','Copy')
%             movefile(outputfile, resultfile );


        end
end
